T1 <-  VGAM::yeo.johnson(Y1, .5, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.5, inverse = T)
pskew <- density(T1, n = 160, from = -4, to = 4)
noskew <- density(T2, n = 160, from = -4, to = 4)
nskew <- density(T3, n = 160, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
pskew
vert <- cbind(pskew$y,noskew$y,nskew$y)
vert
vert <- cbind(seq(-4,4,by = .05),pskew$y,noskew$y,nskew$y)
seq(-4,4,by = .05)
max(pskew$x)
length(pskew$x)
length(seq(-4,4,by = .05))
summary(pskew$x)
summary(seq(-4,4,by = .05))
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
vert
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1.5)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1.5)
T1 <-  VGAM::yeo.johnson(Y1, .5, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.5, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
vert
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .5, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.5, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .25, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.75, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .25, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 3, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .25, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.25, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .25, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.5, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, .90)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .1, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.5, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, 0, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 0, 1)
T1 <-  VGAM::yeo.johnson(Y1, .75, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.25, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, 0, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 0, 1)
T1 <-  VGAM::yeo.johnson(Y1, .50, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.50, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .50, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.50, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
write.csv(vert,"~/desktop/vert.csv")
set.seed(90291)
Y1 <- rnorm(10000000, -2, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 2, 1)
T1 <-  VGAM::yeo.johnson(Y1, .75, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.25, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
write.csv(vert,"~/desktop/vert.csv")
set.seed(90291)
Y1 <- rnorm(10000000, -3, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 3, 1)
T1 <-  VGAM::yeo.johnson(Y1, .50, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.50, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
write.csv(vert,"~/desktop/vert.csv")
set.seed(90291)
Y1 <- rnorm(10000000, 0, 1)
Y2 <- rnorm(10000000, 0, 1)
Y3 <- rnorm(10000000, 0, 1)
T1 <-  VGAM::yeo.johnson(Y1, .50, inverse = T)
T2 <-  VGAM::yeo.johnson(Y2, 1, inverse = T)
T3 <-  VGAM::yeo.johnson(Y3, 1.50, inverse = T)
pskew <- density(T1, n = 161, from = -4, to = 4)
noskew <- density(T2, n = 161, from = -4, to = 4)
nskew <- density(T3, n = 161, from = -4, to = 4)
plot(pskew)
plot(noskew)
plot(nskew)
vert <- cbind(pskew$x,pskew$y,noskew$y,nskew$y)
write.csv(vert,"~/desktop/vert.csv")
install.packages(c("admisc", "bayestestR", "bit", "brew", "broom", "bslib", "Cairo", "callr", "CDM", "class", "cli", "clipr", "clubSandwich", "cluster", "colorspace", "commonmark", "copula", "cpp11", "crayon", "curl", "data.table", "datawizard", "DBI", "desc", "devtools", "digest", "dplyr", "DT", "evaluate", "extrafont", "fansi", "farver", "fastmatrix", "fontawesome", "forcats", "foreign", "furrr", "future", "future.apply", "generics", "gert", "ggplot2", "gh", "gitcreds", "globals", "glue", "gtable", "haven", "hms", "htmltools", "httpuv", "httr", "insight", "isoband", "ISR", "jomo", "jsonlite", "knitr", "lifecycle", "lme4", "magrittr", "markdown", "MASS", "Matrix", "matrixcalc", "MatrixModels", "mclust", "mgcv", "mice", "miceadds", "minqa", "mnormt", "moments", "nlme", "nloptr", "nnet", "openssl", "openxlsx", "ordinal", "parallelly", "parameters", "performance", "pillar", "pkgload", "plyr", "processx", "ps", "psych", "purrr", "r2mlm", "randtoolbox", "RColorBrewer", "Rcpp", "RcppArmadillo", "RcppEigen", "readr", "rlang", "rmarkdown", "rngWELL", "rockchalk", "roxygen2", "rpart", "rprojroot", "rstudioapi", "Rttf2pt1", "rversions", "saemix", "sandwich", "sass", "scales", "scatterplot3d", "shiny", "sirt", "spatial", "stringi", "stringr", "survival", "sys", "sysfonts", "TAM", "testthat", "tibble", "tidyr", "tidyselect", "tinytex", "tzdb", "ucminf", "usethis", "vctrs", "VGAM", "viridisLite", "vroom", "waldo", "withr", "xfun", "XML", "yaml", "zip", "zoo"))
install.packages("ordinal")
install.packages("ordinal")
# example 5.1: curvilinear regression with fiml estimation
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("mathachievement.dat", na.strings = "999")
names(data) <- c("id", "condition", "male", "frlunch", "atrisk", "stanread", "efficacy", "anxiety", "mathpre", "mathpost")
##################################################
# estimate means and center variables
##################################################
# estimate sample statistics in lavaan
model <- "stanread ~ 1; atrisk ~ 1; mathpost ~ 1; anxiety ~ 1; frlunch ~ 1; efficacy ~ 1; mathpre ~ 1;"
descriptives <- inspectSampleCov(model, data, missing = "fiml")
# center lower-order variable
data$anxiety.cgm <- data$anxiety - descriptives$mean["anxiety"]
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread","atrisk","mathpost","anxiety.cgm","frlunch","efficacy")])
# set ranges (nodes) for pseudo-imputations
nodes.atrisk <- c(0,1)
nodes.stanread <- c(20, 80, by = 2)
nodes.mathpost <- seq(30, 90, by = 2)
nodes.anxiety <- seq(-30, 30, by = 2)
nodes.frlunch <- c(0,1)
nodes.efficacy <- seq(1, 6, by = .20)
# model for frlunch predictor: f(frlunch|condition,mathpre)
model.frlunch <- list("model" = "logistic", "formula" = frlunch ~ mathpre, nodes = nodes.frlunch)
# model for efficacy predictor: f(efficacy|frlunch,condition,mathpre)
model.efficacy <- list("model" = "linreg", "formula" = efficacy ~ frlunch + mathpre, nodes = nodes.efficacy)
# model for anxiety predictor: f(anxiety|efficacy,frlunch,condition,mathpre)
model.anxiety <- list("model"="linreg", "formula" = anxiety.cgm ~ efficacy + frlunch + mathpre, nodes = nodes.anxiety)
# model for mathpost outcome: f(mathpost|anxiety,efficacy,frlunch,condition,mathpre)
model.mathpost <- list("model" = "linreg", "formula" = mathpost ~ anxiety.cgm + I(anxiety.cgm^2) + efficacy + frlunch + mathpre, nodes = nodes.mathpost)
# model for atrisk auxiliary variable: f(atrisk|mathpost,efficacy,frlunch,condition,mathpre)
model.atrisk <- list("model" = "logistic", "formula" = atrisk ~ mathpost + anxiety.cgm + efficacy + frlunch + mathpre, nodes = nodes.atrisk)
# model for stanread auxiliary variable: f(stanread|atrisk,mathpost,efficacy,frlunch,condition,mathpre)
model.stanread <- list("model" = "linreg", "formula" = stanread ~ atrisk + mathpost + anxiety.cgm + efficacy + frlunch + mathpre, nodes = nodes.stanread)
# combine predictor models into a list
predictor.models <- list(frlunch = model.frlunch, efficacy = model.efficacy, anxiety = model.anxiety, mathpost = model.mathpost, atrisk = model.atrisk)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread, ind = predictor.models)
summary(fit)
# example 5.1: curvilinear regression with fiml estimation
library(fdir)
library(lavaan)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("mathachievement.dat", na.strings = "999")
names(data) <- c("id", "condition", "male", "frlunch", "atrisk", "stanread", "efficacy", "anxiety", "mathpre", "mathpost")
##################################################
# estimate means and center variables
##################################################
# estimate sample statistics in lavaan
model <- "stanread ~ 1; atrisk ~ 1; mathpost ~ 1; anxiety ~ 1; frlunch ~ 1; efficacy ~ 1; mathpre ~ 1;"
descriptives <- inspectSampleCov(model, data, missing = "fiml")
# center lower-order variable
data$anxiety.cgm <- data$anxiety - descriptives$mean["anxiety"]
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread","atrisk","mathpost","anxiety.cgm","frlunch","efficacy")])
# set ranges (nodes) for pseudo-imputations
nodes.atrisk <- c(0,1)
nodes.stanread <- c(20, 80, by = 2)
nodes.mathpost <- seq(30, 90, by = 2)
nodes.anxiety <- seq(-30, 30, by = 2)
nodes.frlunch <- c(0,1)
nodes.efficacy <- seq(1, 6, by = .20)
# model for frlunch predictor: f(frlunch|condition,mathpre)
model.frlunch <- list("model" = "logistic", "formula" = frlunch ~ mathpre, nodes = nodes.frlunch)
# model for efficacy predictor: f(efficacy|frlunch,condition,mathpre)
model.efficacy <- list("model" = "oprobit", "formula" = efficacy ~ frlunch + mathpre, nodes = nodes.efficacy)
# model for anxiety predictor: f(anxiety|efficacy,frlunch,condition,mathpre)
model.anxiety <- list("model"="linreg", "formula" = anxiety.cgm ~ efficacy + frlunch + mathpre, nodes = nodes.anxiety)
# model for mathpost outcome: f(mathpost|anxiety,efficacy,frlunch,condition,mathpre)
model.mathpost <- list("model" = "linreg", "formula" = mathpost ~ anxiety.cgm + I(anxiety.cgm^2) + efficacy + frlunch + mathpre, nodes = nodes.mathpost)
# model for atrisk auxiliary variable: f(atrisk|mathpost,efficacy,frlunch,condition,mathpre)
model.atrisk <- list("model" = "logistic", "formula" = atrisk ~ mathpost + anxiety.cgm + efficacy + frlunch + mathpre, nodes = nodes.atrisk)
# model for stanread auxiliary variable: f(stanread|atrisk,mathpost,efficacy,frlunch,condition,mathpre)
model.stanread <- list("model" = "linreg", "formula" = stanread ~ atrisk + mathpost + anxiety.cgm + efficacy + frlunch + mathpre, nodes = nodes.stanread)
# combine predictor models into a list
predictor.models <- list(frlunch = model.frlunch, efficacy = model.efficacy, anxiety = model.anxiety, mathpost = model.mathpost, atrisk = model.atrisk)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread, ind = predictor.models)
summary(fit)
# example 2.1:logistic regression with fiml estimation
library(fdir)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","atriskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","read9grp","stanread7",
"math1","math2","math3","math9","math9grp","stanmath7")
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9grp","read1","lrnprob1","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.hispanic <- c(0,1)
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9grp <- c(0,1)
nodes.read1 <- seq(30, 160, by = 2)
nodes.lrnprob1 <- seq(20, 100, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor: f(behsymp1)
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor: f(lrnprob1|behsymp1)
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1 ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor: f(read1|lrnprob1,behsymp1)
model.read1 <- list( "model" = "linreg", "formula" = read1 ~ lrnprob1 + behsymp1, nodes = nodes.read1)
# model for read9grp outcome: f(read9grp|read1,lrnprob1,behsymp1)
model.read9grp <- list( "model" = "logistic", "formula" = read9grp ~ read1 + lrnprob1 + behsymp1, nodes = nodes.read9grp)
# model for hispanic auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.hispanic <- list("model" = "logistic", "formula" = hispanic ~ read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.hispanic)
# model for read2 auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ hispanic + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable f(stanread7|read2,read9grp,read1,lrnprob1,behsymp1)
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + hispanic + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(read2 = model.read2, hispanic = model.hispanic, read9grp = model.read9grp, behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# example 2.1:logistic regression with fiml estimation
library(fdir)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","read9grp","stanread7",
"math1","math2","math3","math9","math9grp","stanmath7")
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9grp","read1","lrnprob1","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9grp <- c(0,1)
nodes.read1 <- seq(30, 160, by = 2)
nodes.lrnprob1 <- seq(20, 100, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor: f(behsymp1)
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor: f(lrnprob1|behsymp1)
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1 ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor: f(read1|lrnprob1,behsymp1)
model.read1 <- list( "model" = "linreg", "formula" = read1 ~ lrnprob1 + behsymp1, nodes = nodes.read1)
# model for read9grp outcome: f(read9grp|read1,lrnprob1,behsymp1)
model.read9grp <- list( "model" = "logistic", "formula" = read9grp ~ read1 + lrnprob1 + behsymp1, nodes = nodes.read9grp)
# model for read2 auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable f(stanread7|read2,read9grp,read1,lrnprob1,behsymp1)
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(read2 = model.read2, read9grp = model.read9grp, behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# example 2.1:logistic regression with fiml estimation
library(fdir)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","read9grp","stanread7",
"math1","math2","math3","math9","math9grp","stanmath7")
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read9grp","read1","lrnprob1","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9grp <- c(0,1)
nodes.read1 <- seq(30, 160, by = 2)
nodes.lrnprob1 <- seq(20, 100, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor: f(behsymp1)
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor: f(lrnprob1|behsymp1)
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1 ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor: f(read1|lrnprob1,behsymp1)
model.read1 <- list( "model" = "linreg", "formula" = read1 ~ lrnprob1 + behsymp1, nodes = nodes.read1)
# model for read9grp outcome: f(read9grp|read1,lrnprob1,behsymp1)
model.read9grp <- list( "model" = "logistic", "formula" = read9grp ~ read1 + lrnprob1 + behsymp1, nodes = nodes.read9grp)
# model for read2 auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable f(stanread7|read2,read9grp,read1,lrnprob1,behsymp1)
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1, read9grp = model.read9grp, read2 = model.read2)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
# FIML Logistic Regression with a Binary Outcome
library(fdir)
library(mdmb)
##################################################
# read data
##################################################
# set working directory
set()
# read raw data from working directory
data <- read.table("behaviorachievement.dat", na.strings = "999")
names(data) <- c("id","male","hispanic","riskgrp","atrisk","behsymp1","lrnprob1",
"read1","read2","read3","read9","read9grp","stanread7",
"math1","math2","math3","math9","math9grp","stanmath7")
##################################################
# analysis
##################################################
# summarize incomplete predictors to determine ranges for pseudo-imputations
summary(data[,c("stanread7","read2","read1","lrnprob1","behsymp1")])
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9grp <- c(0,1)
nodes.read1 <- seq(30, 160, by = 2)
nodes.lrnprob1 <- seq(20, 100, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor:
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1 ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor
model.read1 <- list( "model" = "linreg", "formula" = read1 ~ lrnprob1 + behsymp1, nodes = nodes.read1)
# model for read9grp outcome
model.read9grp <- list( "model" = "logistic", "formula" = read9grp ~ read1 + lrnprob1 + behsymp1, nodes = nodes.read9grp)
# focal model for read2 auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1, read9grp = model.read9grp, read2 = model.read2)
# estimate factored regression model w mdmb
fit <- frm_em(dat = data, dep = model.stanread7, ind = predictor.models)
summary(fit)
library(mdmb)
data_url <- "https://raw.githubusercontent.com/craigenders/ies-toolkit/main/Data/behaviorachievement.rda"
load(gzcon(url(data_url, open = "rb")))
# summaries to determine ranges for pseudo-imputations
summary(behaviorachievement)
# summaries to determine ranges for pseudo-imputations
summary(behaviorachievement)
# set ranges (nodes) for pseudo-imputations
nodes.stanread7 <- seq(80, 420, by = 5)
nodes.read2 <- seq(10, 160, by = 2)
nodes.read9grp <- c(0,1)
nodes.read1 <- seq(30, 160, by = 2)
nodes.lrnprob1 <- seq(20, 100, by = 2)
nodes.behsymp1 <- seq(10, 100, by = 2)
# model for behsymp1 predictor
model.behsymp1 <- list( "model" = "linreg", "formula" = behsymp1 ~ 1, nodes = nodes.behsymp1)
# model for lrnprob1 predictor:
model.lrnprob1 <- list( "model" = "linreg", "formula" = lrnprob1 ~ behsymp1, nodes = nodes.lrnprob1)
# model for read1 predictor
model.read1 <- list( "model" = "linreg", "formula" = read1 ~ lrnprob1 + behsymp1, nodes = nodes.read1)
# model for read9grp outcome
model.read9grp <- list( "model" = "logistic", "formula" = read9grp ~ read1 + lrnprob1 + behsymp1, nodes = nodes.read9grp)
# focal model for read2 auxiliary variable f(read2|read9grp,read1,lrnprob1,behsymp1)
model.read2 <- list("model" = "linreg", "formula" = read2 ~ read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.read2)
# model for stanread7 auxiliary variable
model.stanread7 <- list("model" = "linreg", "formula" = stanread7 ~ read2 + read9grp + read1 + lrnprob1 + behsymp1, nodes = nodes.stanread7)
# combine predictor models into a list
predictor.models <- list(behsymp1 = model.behsymp1, lrnprob1 = model.lrnprob1, read1 = model.read1, read9grp = model.read9grp, read2 = model.read2)
# estimate factored regression model w mdmb
fit <- frm_em(dat = behaviorachievement, dep = model.stanread7, ind = predictor.models)
summary(fit)
summary(fit)
